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Abstract 

Algorithms are presented for the tanh- and sech-methods, which lead to 
closed-form solutions of nonlinear ordinary and partial differential equa- 
tions (ODEs and PDEs). New algorithms are given to find exact polyno- 
mial solutions of ODEs and PDEs in terms of Jacobi's elliptic functions. 

For systems with parameters, the algorithms determine the conditions 
on the parameters so that the differential equations admit polynomial 
solutions in tanh, sech, combinations thereof, Jacobi's sn or cn functions. 
Examples illustrate key steps of the algorithms. 

The new algorithms are implemented in Mathematica. The package 
PDESpecialSolutions.m can be used to automatically compute new spe- 
cial solutions of nonlinear PDEs. Use of the package, implementation 
issues, scope, limitations, and future extensions of the software are ad- 
dressed. 

A survey is given of related algorithms and symbolic software to com- 
pute exact solutions of nonlinear differential equations. 
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1. Introduction 

The appearance of solitary wave solutions in nature is quite common. Bell shaped 
sech-solutions and kink shaped tanh-solutions model wave phenomena in fluids, 
plasmas, elastic media, electrical circuits, optical fibers, chemical reactions, bio- 
genetics, etc. The travelling wave solutions of the Korteweg-de Vries (KdV) and 
Boussinesq equations, which describe water waves, are famous examples. 

Apart from their physical relevance, the knowledge of closed-form solutions of 
nonlinear ordinary and partial differential equations (ODEs and PDEs) facili- 
tates the testing of numerical solvers, and aids in the stability analysis. Indeed, 
the exact solutions given in this paper correspond to homoclinic and heteroclinic 
orbits in phase space, which are the separatrices of stable and unstable regions. 

Travelling wave solutions of many nonlinear ODEs and PDEs from soliton 
theory (and beyond) can often be expressed as polynomials of the hyperbolic 
tangent and secant functions. An explanation is given in, for example, Hereman 
and Takaoka (1990). The existence of solitary wave solutions of evolution equa- 
tions is addressed in Kichenassamy and Olver (1993). The tanh-method provides 
a straightforward algorithm to compute such particular solutions for a large class 
of nonlinear PDEs. Consult Malfliet (1992, 2003), Malfliet and Hereman (1996), 
and Das and Sarma (1999) for a multitude of references to tanh-based techniques 
and applications. 

The tanh-method for, say, a single PDE in u(x, t) works as follows: In a travel- 
ling frame of reference, £ = c\x + c^t + A, one transforms the PDE into an ODE 
in the new independent variable T = tanh£. Since the derivative of tanh is poly- 
nomial in tanh, i.e., T' — 1 — T 2 , all derivatives of T are polynomials of T. Via a 
chain rule, the polynomial PDE in u(x,t) is transformed into an ODE in U(T), 
which has polynomial coefficients in T. One then seeks polynomial solutions of 
the ODE, thus generating a subset of the set of all solutions. 

Along the path, one encounters ODEs which are nonlinear, higher-order ver- 
sions of the ultraspherical differential equation, 

(1 - x 2 )y"{x) - {2a + l)xy'{x) + n(n + 2a)y{x) = 0, (1) 

with integer n > and a real, whose solutions are the Gegenbauer polynomials. 
Eq. (JTJ) includes the Legendre equation (a = ~), satisfied by the Legendre poly- 
nomials, and the ODEs for Chebeyshev polynomials of type I (a = 0) and type 
II (a = 1). Likewise, the associated Legendre equation, 

(1 - x 2 ) 2 y"{x) - 2x[x 2 - l)y'(x) + [n{n + 1)(1 - x 2 ) - m 2 ]y(x) = 0, (2) 

with m and n non-negative integers, appears in solving the Sturm-Liouville prob- 
lem for the KdV with a sech-square potential (see Drazin and Johnson, 1989). 

The appeal and success of the tanh-method lies in the fact that one circum- 
vents integration to get explicit solutions. Variants of the method appear in 
mathematical physics, plasma physics, and fluid dynamics. For early references 



D. Baldwin et al.: Symbolic computation of Tanh, Sech, Cn, and Sn Solutions 3 

see e.g. Malfliet, 1992; Yang, 1994; and Das and Sarma, 1999. Recently, the 
tanh-methods have been applied to many nonlinear PDEs in multiple indepen- 
dent variables (see Fan, 2002abc, 2003abc; Fan and Hon, 2002, 2003ab; Gao and 
Tian, 2001; Li and Liu, 2002; Yao and Li, 2002ab). 

In this paper we present three flavors of tanh- and sech-methods as they apply 
to nonlinear polynomial systems of ODEs and PDEs. Based on the strategy of 
the tanh-method, we also present algorithms to compute polynomial solutions 
in terms of the Jacobi sn and cn functions. Applied to the KdV equation, the 
so-called cnoidal solution (Drazin and Johnson, 1989) is obtained. For Duffing's 
equation (Lawden, 1989), we recover known sn and cn-solutions which model 
vibrations of a nonlinear spring. Sn- and cn-methods are quite effective for sym- 
bolically solving nonlinear PDEs as shown in Fu et al. (2001), Parkes et al. 
(2002), Liu and Li (2002ab), Fan and Zhang (2002), Fan (2003abc), Chen and 
Zhang (2003ab), and Yan (2003). 

We also present our package, PDESpecialSolutions .m (Baldwin et al, 2001) 
in Mathematica, which implements the five methods. Without intervention by 
the user, our software computes travelling wave solutions as polynomials in either 
T = tanh£, S = sech£, combinations thereof, CN = cn(£; m), or SN = sn(£; m) 
with £ = C\X + c%y + c%z + . . . + c n t + A = J2f=o c j x j + A. The coefficients of the 
spatial coordinates are the components of the wave vector; the time coefficient 
is the angular frequency of the wave. The wave travels in the direction of the 
wave vector; its plane wave front is perpendicular to that wave vector. A is the 
constant phase. For systems of ODEs or PDEs with constant parameters, the 
software automatically determines the conditions on the parameters so that the 
equations might admit polynomial solutions in tanh, sech, both, sn or cn . 

Parkes and Duffy (1996) mention the difficulty of using the tanh-method by 
hand for anything but simple PDEs. Therefore, they automated to some degree 
the tanh-method using Mathematica. Their code ATFM carries out some (but 
not all) steps of the method. Parkes et al. (1998) also considered solutions to 
(odd-order generalized KdV) equations in even powers of sech . The code ATFM 
does not cover solutions involving odd powers of sech . Recently, Parkes et al. 
(2002) extended their methods to cover the Jacobi elliptic functions. Abbott et al. 
(2002) produced the function SeriesSn to partially automate the elliptic function 
method. Li and Liu (2002) designed the Maple package RATH to automate the 
tanh-method. In Liu and Li (2002a) they announce their Maple code AJFM for 
the Jacobi elliptic function method. In Section I5~2l we review the codes ATFM, 
RATH, AJFM, and SeriesSn and compare them with PDESpecialSolutions. m. 

The paper is organized as follows: In Sections 121 and El we give the main steps 
of the algorithms for computing tanh- and sech-solutions of nonlinear polynomial 
PDEs. We restrict ourselves to polynomial solutions in either tanh or sech . The 
Boussinesq equation and Hirota-Satsuma system of coupled KdV equations illus- 
trate the steps. For references to both equations see e.g. Ablowitz and Clarkson 
(1991). In Section0]we consider a broader class of polynomial solutions involving 
both tanh and sech . The tanh-sech algorithm is used to solve a system of PDEs 
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due to Gao and Tian (2001). In Section El we show how modifying the chain 
rule allows us to find polynomial solutions in cn and sn. The KdV equation is 
used to illustrates the steps. In Section |H1 we give details of the algorithms to 
compute the highest-degree of the polynomials, to analyze and solve nonlinear 
algebraic systems with parameters, and to numerically and symbolically test so- 
lutions. The coupled KdV equations illustrate the subtleties of these algorithms. 
In Section [7| we present exact solutions for several nonlinear ODEs and PDEs. 
In Section |S1 we address other perspectives and extensions of the algorithms, and 
review related software packages. We discuss the results and draw some conclu- 
sions in Section El The use of the package PDESpecialSolutions.m is shown in 
Appendix A. 



2. Algorithm to compute tanh-solutions for nonlinear PDEs 

In this section we outline the tanh-method (Malfliet and Hereman, 1996) for 
the computation of closed- form tanh-solutions for nonlinear PDEs (and ODEs). 
Each of the five main steps of our algorithm is illustrated for the Boussinesq 
equation. Details of steps T2, T4 and T5 are postponed to Sectional 
Given is a system of polynomial PDEs with constant coefficients, 

A(u(x), u'(x), u"(x), ■ • • , u«(x), • ■ • , uW(x)) = 0, (3) 

where the dependent variable u has M components Ui, the independent variable 
x has N components Xj, and u( fc )(x) denotes the collection of mixed derivative 
terms of order k. Lower-case Greek letters will denote parameters in (J3J). 

For notational simplicity, in Section [7| we will use dependent variables u,v,w, 
etc. and independent variables x, y, z, and t. 
Example: The classical Boussinesq equation, 

u t t - u xx + 3uu xx + 3u 2 x + <\ti,,,, = 0, (4) 

with real parameter a, was proposed by Boussinesq to describe surface wa- 
ter waves whose horizontal scale is much larger than the depth of the water 
(Ablowitz and Clarkson, 1991). Variants of (j3J) were recently solved by Fan and 
Hon (2003a). 

While one could apply the tanh-method directly to (UJ), we recast it as a first 
order system in time to show the method for a simple system of PDEs. So, 

Ui, X2 + u 2, Xl = 0, 

U2,x 2 + v>i,xi ~ 3 u iUi, xl ~ au 1)3xi = 0, 
where x\ = x,X2 = t, U\(xx, x-z) = u(x, t), and 1*2(^1, £2) = u t (x, t). We use 

def d k Uj 

U iMj — Q x k ! 

through out this paper. 



(5) 



U; 



def 



r+s 



U: 



,pxj rx k SX( 



dx V jdx r k dx\ ' 



etc. 



(6) 
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Step Tl: Transform the PDE into a nonlinear ODE 

We seek solutions in the travelling frame of reference, 

N 

where Cj and A are constant. 

The tanh-method seeks polynomial solutions expressible in the hyperbolic 
tangent, T = tanh£. Based on the identity cosh 2 £ — sinh 2 £ = 1 one computes 

tanh'f = sech 2 £ = 1 — tanh 2 £ , 
tanh"f = -2tanh£ + 2tanh 3 £, etc. 

Therefore, the first and, consequently, all higher-order derivatives are polynomi- 
als in T. Since T" = 1 — T 2 , repeatedly applying the chain rule, 

<h_ chdTd£_ _ 2 dm_ 

ax, dTd^dxj A Ur K) 

transforms the system of PDEs into a coupled system of nonlinear ODEs, 

A(T, U(T), U'(T), U"(T), . . . , U (m) (T)) = 0, (10) 

with U(T) = u(x). Each component of A is a nonlinear ODE with polynomial 
coefficients in T. 
Example: Substituting 

u hX] =c 3 (l-T 2 )Ui 

u ij2xj = c 2 (l - T 2 ) [(1 - T 2 )f/;]' = c 2 (l - T 2 )\-2TU[ + (1 - T 2 )U"], 
u ijSxj = cf(l - T 2 ) [-2T(1 - T 2 )U- + (1 - T 2 ) 2 U'^ ' ' ' ' 

= c 3 (l - T 2 ) [-2(1 - 3T 2 )f/; - 6T(1 - T 2 )U" + (1 - T 2 ) 2 ^"] , 

into (0), and cancelling common (1 — T 2 ) factors, yields 
c 2 U[ + Cl U' 2 = 0, 

c^+^UiS^UtUi+acl [2(1-3T 2 )U[+6T(1-T 2 )U^-(1-T 2 ) 2 U^] =0, (12) 

where U\(T) = u%(xi,x 2 ) and U 2 (T) = u 2 (xi,x 2 ). 

Step T2: Determine the degree of the polynomial solutions 

Seeking polynomial solutions of the form 

Mi 

cMT) = 5>,T J , (13) 

j=0 

we must determine the leading exponents Mi before the a^- can be computed. 



D. Baldwin et al.: Symbolic computation of Tanh, Sech, Cn, and Sn Solutions 6 



We assume that Mj > 1 to avoid trivial solutions. Substituting Ui into (|TU|) . the 
coefficients of every power of T in every equation must vanish. In particular, the 
highest degree terms must vanish. Since the highest degree terms depend only on 
T Mi in (0, it suffices to substitute Ui(T) = T M% into the left hand side of ([11])). 
In the resulting polynomial system P(T), equating every two possible highest 
exponents in every component Pj gives a linear system for the Mj. That linear 
system is then solved. 

If one or more exponents Mj remain undetermined, assign an integer value to 
the free Mj so that every equation in (jlUj) has at least two different terms with 
equal highest exponents. Carry each solution to step T3. 

Example: For the Boussinesq system, substituting U\{T) = T Ml and U 2 (T) = 
T Mz into (j!2|) . and equating the highest exponents of T for each equation, gives 

Mi - 1 = M 2 - 1, 2Mi - 1 = Mi + 1. (14) 

Then, M x = M 2 = 2, and 

U X {T) = a w + a n T + a 12 T 2 , U 2 {T) = a 20 + a 21 T + a 22 T 2 . (15) 

Step T3: Derive the algebraic system for the coefficients 

To generate the system for the unknown coefficients and wave parameters 
Cj, substitute (fT3Jl into (jTUj) and set the coefficients of T % to zero. The resulting 
nonlinear algebraic system for the unknowns a^- is parameterized by the Cj, and 
the external parameters (in lower-case Greek letters) of system (J3j), if any. 
Example: Continuing with the Boussinesq system, after substituting ({To} into 
(|12p. and collecting the terms of like degree in T, we get (in order of complexity): 



a 21 Ci + an c 2 = 0, 
a 22 ci + ai2 c 2 = 0, 
an Ci (3ai2 + 2a cf) = 0, 
ai2Ci (ai2 + 4a cf) = 0, 
an Ci — 3aio an c\ + 2aan c\ + 021 C2 = 0, 
-3a n ci + 2 012 c x — 6ai a 12 c\ + 16a a i2 c\ + 2a 2 2 c 2 = 0, 



(16) 



with unknowns ai , an, a i2 , a 20 , a 2 i, a 2 2, and parameters Ci, c 2 , and a. 
Step T4: Solve the nonlinear parameterized algebraic system 

The most difficult step is solving the nonlinear algebraic system. To do so, we 
designed a customized, yet powerful, nonlinear solver (see Section l6~21 for details). 

The nonlinear algebraic system is solved with the following assumptions: 
(i) all parameters, a, [3, etc., in (jSJ) are strictly positive. Vanishing parameters 
may change the exponents Mj in step T2. To compute solutions corresponding 
to negative parameters, reverse the signs of the parameters in the PDE. For 
example, replace a by —a in (J3J). 
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(ii) the coefficients of the highest power terms (a,iMi, i = 1, • • • , M) in (fT^J) are 
all nonzero (for consistency with step T2). 

(iii) all Cj are nonzero (demanded by the physical nature of the solutions). 
Example: Assuming c\, c 2 , «i2, 022, and a nonzero, the solution of is 



In this case, there are no conditions on the parameters ci, c 2 and a. 
Step T5: Build and test the solitary wave solutions 

Substitute the solutions obtained in step T4 into (fT3*j) and reverse step Tl to 
obtain the explicit solutions in the original variables. It is prudent to test the 
solutions by substituting them into ©. For details about testing see Section 

Example: Inserting ()17|) into (j!5|) . and replacing T = tanh(cix + c 2 t + A), the 
closed form solution for © (or (jlj)) is 



where 020, ci, C2, a and A are arbitrary. Steps T1-T5 must be repeated if one or 
more of the external parameters (lower-case Greeks) are set to zero. 

3. Algorithm to compute sech-solutions for nonlinear PDEs 

In this section we restrict ourselves to polynomial solutions of (jHJ) in sech . Poly- 
nomial solutions involving both sech and tanh are dealt with in Section 0] Details 
of the algorithms for steps S2, S4 and S5 are in Section El 

Using tanh 2 £ + sech 2 £ = 1, solution (jTSjl of © can be expressed as 

Ui(x,t) = (c 2 — c 2 — 4ac^)/(3c 2 ) + 4ac 2 sech 2 (cix + c 2 t + A), 

u 2 (x, t) = a 20 + 4aciC2 — 4aciC2 sech 2 (cix + c^t + A). ^ 

Obviously, any even order solution in tanh can be written in even orders of 
sech . Some PDEs however have polynomial solutions of odd-order in sech . For 
example, the modified KdV equation (Ablowitz and Clarkson, 1991), 



a io = ( c i ~ c 2 + 8 ac i)/(3c 1 ), an = 0, (X12 = —Aac^ 
a2o = arbitrary, 021 = 0, a 2 2 = 4aciC2. 



(17) 



u(x, t) = Ui(x,t) = (c 2 — c\ + 8ac^)/(3c 2 ) — 4ac 2 tanh 2 (cix + c^t + A) 
u 2 (x, t) = —J ui ;t (x, t)dx = a 2 o + 4aciC2 tanh 2 (cix + c 2 t + A), 



(18) 



u t + au 2 u x + u. 



'XXX 



(20) 



has the solution 



u 



(x, t) = ±Ci \/6/asech(cia; — c\t + A) 



(21) 



which cannot be found using the tanh-method. 

Example: The five main steps of the sech-algorithm are illustrated with the 
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Hirota-Satsuma system of coupled KdV equations (Ablowitz and Clarkson, 1991), 

u,-o(6».+«_)+!!**. = 0, 
v t + 3uv x + v xxx = 0, 

with real parameters a,/3. Sech-type solutions were reported in Hereman (1991) 
and Fan and Hon (2002). Variants and generalizations of (|2*2*|) were solved in 
Chen and Zhang (2003a) and Yan (2003). 

Letting ui(xi,x 2 ) = u(x,t) and u 2 (xi,x 2 ) = v(x,t), Eq. (|2*2*Jl is then 

u 1>X2 - a(6wiUi )Xl + ux i3xi ) + 2[3u 2 u 2 , Xl =0, 

U2,x 2 +3MiM2,X! +^2,3^ =0. 

Step SI: Transform the PDE into a nonlinear ODE 

Adhering to the travelling frame of reference (J7J), and using tanh 2 £ + sech 2 £ = 1, 



sech'£ = — sech£ tanh£ = — sech£y 1 — sech 2 £. (24) 
Setting S = sech £ and repeatedly applying the chain rule, 

* = *L§ <* = -c jS VT^% (25) 
dxj dS d£ GXj dS 

(JHJ) is transformed into a system of nonlinear ODEs of the form: 



T(S, U(S), U'(S), ...) + VT^S* U(S, V(S), V'(S), . . .) = 0, (26) 

where U(5) = u(x), and all components of T and II are ODEs with polynomial 
coefficients in S. If either T or II is identically 0, then 

A(S,U(S),U / (S),...) = 0, (27) 

where A is either T or II, whichever is nonzero. For this to occur, the order of 
all terms in any equation in © must be even or odd (as is the case in ([23)1 ). 

Any term in (J3J) for which the total number of derivatives is even contributes 
to the first term in ()2fij) : while any term of odd order contributes to the second 
term. Section 0] deals with any case for which neither T or II is identically 0. 
Example: Substituting 



iSVT=S*UL 



u 



l,XjX k — 



CiCkSVT^^s* SVT^PU' 



= c jCk S[(l - 2S 2 )U[ + S(l - S 2 )UH (28) 

»„,.„.,-, = -CiwsVT^s* [s(i - 2s 2 )u[ + s(i - s 2 )u';]' 



—c 



•.j c kCl S VT^P [( 1 - QS 2 ) U[ + 3S ( 1 - 2S 2 ) U'l + S 2 ( 1 - S 2 ) U?'] 
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into (J2HJ), and cancelling the common Sy/l — S 2 factors yields 

c 2 U[-6a Cl U x U[-ac\ [(1-6S 2 ) U[+3S( 1-2S 2 ) U'(+S 2 ( 1-S 2 ) U"']+2^ Cl U 2 U 2 ^ 0, 
c 2 c^+3 Cl cW2+c? [(l-QS 2 )U 2 +3S(l-2S 2 )U'^+S 2 (l-S 2 )U' 2 "] = 0, (29) 

with Ui(T) = Ui(xi,x 2 ) and U 2 (T) = u 2 (xi,x 2 ). Note that (J2T?)) matches (|27j) 
with A = II, since T = 0. 

Step S2: Determine the degree of the polynomial solutions 

We seek polynomial solutions of the form, 

Mi 

U i (S) = J2a ij S j . (30) 

j=0 

To determine the Mj exponents, substitute Ui(S) = S Mi into the left hand side 
of (|27j) and proceed as in step T2. Continue with step S3 for each solution of 
Mj. If some of the Mj exponents are undetermined, try all legitimate values for 
the free Mj. See Section loTTl for more details. 

Example: For (J23J), substituting U X {S) = S M \U 2 {S) = S m into flU and 
equating the highest exponents in the second equation yields Mi + M 2 — 1 = 
1 + M 2 , or Mi = 2. The maximal exponents coming from the first equation are 
2Mi - 1 (from the V X V' X term), Mi + 1 (from U'{' ), and 2M 2 - 1 (from U 2 U' 2 ). 
Using Mi = 2, two cases emerge: (i) the third exponent is less than the first two 
(equal) exponents, i.e., 2M 2 — 1 < 3, so M 2 = 1, or (ii) all three exponents are 
equal, in which case M 2 = 2. For the case Mi = 2 and M 2 = 1, 

U(S) = a w + a lx S + a 12 S 2 , V(S) = a 20 + a 21 S, (31) 

and, for the case Mi = M 2 = 2, 

U(S) = a w + auS + ai 2 S 2 , V(S) = a 20 + a 2 iS + a 22 S 2 . (32) 

Step S3: Derive the algebraic system for the coefficients 

Follow the strategy in step T3. 

Example: After substituting (|31|) into (|29|). cancelling common numerical fac- 
tors, and organizing the equations (according to complexity) one obtains 

a n a 2 ic± = 0, 
aaiiCi(3ai2 — cf) = 0, 
aai 2 ci(ai 2 - 2c\) = 0, 
a 21 Cx{a 12 - 2c\) = 0, (33) 
a 2 i(3ai Ci + c\ + c 2 ) = 0, 
6aa 10 anCi — 2/3a 20 a 2 iCi + aa u cl — a\ 2 c 2 = 0, 
Saa^ci + 6aaioai 2 ci — /9a 2 i c i + 4aai 2 Ci — ai 2 c 2 = 0. 
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Similarly, after substitution of (jH2j) into (J2HJ), one gets 

a 22 ci(ai2 - 4c 2 ) = 0, 
a 2 i(3a 10 ci + cl + c 2 ) = 0, 
ci(a 12 a 2 i + 2ana 2 2 — 2a 2 \c\) = 0, 



Ci(3aanai2 — (3a 2 \a 22 — crane 2 ) = 0, (34) 

Ci(3aa 2 2 — /^a 2 ^ — 6aai 2 c 2 ) = 0, 

6aaioanCi — 2/3a 2 oa 2 iCi + aauCj — anC2 = 0, 

3a n a 21 ci + 6a w a 22 ci + Sa 22 c\ + 2a 22 c 2 = 0, 

3aa n Ci + 6aa 2 a 12 ci — /3a 2 i c i — 2/3a 20 a 22 Ci + 4aai 2 c5 — a 12 c 2 = 0. 

Step S4: Solve the nonlinear parameterized algebraic system 

Similar strategy as in step T4. 

Example: For a, f3,Ci,c 2 ,a 12 and a 2 \ all nonzero, the solution of ()33|) is 



oio = ~{cl + c 2 )/(3c x ), an = 0, a i2 = 2c 2 , 



a 20 = 0, a 2i = ±y / (4ac^ - 2(1 + 2a)c l c 2 )/f3. 
For a, /?, ci, c 2 , a 12 and a 22 nonzero, the solution of (J5%)l is 
a 10 = -(4ci + c 2 )/ (3ci), an = 0, a i2 = 4c 2 , 

a 20 = ±(4ac^ + (1 + 2a)c 2 )/(ci\/6a/3), a 21 = 0, a 22 = T2c\ \/6a//3. 



(35) 



(36) 



Step S5: Build and test the solitary wave solutions 

Substitute the result of step S4 into (|5U|l and reverse step SI. Test the solutions. 
Example: The solitary wave solutions of ()23|) are 



u (x, t) = -(of + c 2 ) /(3ci) + 2c\ sech 2 (cix + c 2 t + A), 

(37) 



v(x, t) = ±J [Aacf - 2(1 + 2a)dc 2 ]/(3 sech(cia; + c 2 t + A), 
and 

u(x, t) = -(4c? + c 2 ) /(3ci) + 4c 2 sech 2 (cix + c 2 t + A), 

w(x,t) = ±(4ac? + (l + 2a)c 2 )/(ci V / 6a^) =f 2c\ a/ 6a / (3 sech 2 (cia; +c 2 t + A) . ^ 38 ^ 

In both cases ci, c 2 , a, /?, and A are arbitrary. These solutions contain the solu- 
tions reported in Hereman (1991). 

Steps S1-S5 must be repeated if any of the parameters in Q are set to zero. 



D. Baldwin et al: Symbolic computation of Tanh, Sech, Cn, and Sn Solutions 11 

4. Algorithm for mixed tanh-sech solutions for PDEs 

The five main steps of our algorithm to compute mixed tanh-sech solutions for 
(JSJ) are presented below. Here we seek particular solutions of ()26|) when r ^ 
and n ^ 0. On could apply the method of Section El to in 'squared' form 
T 2 (S, U(S), U'(S), . . .) - (1 - S 2 ) U 2 (S, U(5), U'(S), . . .) = 0. For anyth ing but 
simple cases, the computations are unwieldy. Alternatively, since T = Vl — S 2 , 
Eq. may admit solutions of the form 

Mi Ni 

u l (s) = J2H~ a ^ S3Tk - ( 39 ) 

However, can always be rearranged such that 

Mi Ni Ah JVj 

Ui{S) = a l3 S j + T hjSi = a H SJ + v 7 ! 3 ^ b l3 S j . (40) 

The polynomial solutions in S from Section El are special cases of this broader 
class. Remarkably, (j2ZI) where ^/l — S 2 is not explicitly present also admits so- 
lutions of the form ()4U|). See Section ITol for an example. 

Computing solutions of type (|30p with the tanh-sech method is inefficient and 
costly, as the following example and the examples in Sections 17.51 and 17.61 show. 
Example: We illustrate this algorithm with the system (Gao and Tian, 2001): 

u t — u x — 2v = 0, 

v t + 2uw = 0, (41) 
w t + 2uv = 0. 

Step ST1: Transform the PDE into a nonlinear ODE 

Same as step SI. 

Example: Use (|23|) to transform (j4"Tj) into 

(ci - c 2 )SVl - S 2 U[ - 2U 2 = 0, 

c 2 SVl - S 2 U' 2 - 2UxU 3 = 0, (42) 
c 2 SVl - S 2 U' 3 - 2UxU 2 = 0. 

with Ui(S) = Ui(xi, x 2 ), i — 1,2, 3. 

Step ST2: Determine the degree of the polynomial solutions 

Seeking solutions of form ()40|) . we must first determine the leading Mj and A^j 
exponents. Substituting Ui(S) = a ia + a iM .S Mi + y/1 — S 2 (b i0 + b iNi S Ni ) into 
the left hand side of (J2SJ), we get an expression of the form 

P(S) + VT^S 2 Q(S), (43) 



where P and Q are polynomials in S. 
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Consider separately the possible balances of highest exponents in all Pi and Q{. 
Then solve the resulting linear system(s) for the unknowns Mj and iVj. Continue 
with each solution in step ST3. 

In contrast to step S2, we no longer assume Mi > l,iVj > 1. Even with 
some Mi or iVj zero, non-constant solutions Ui(S) often arise. In most examples, 
however, the sets of balance equations for Mj and iVj are too large or the cor- 
responding linear systems are under-determined (i.e., several leading exponents 
remain arbitrary). To circumvent the problem, we set all Mj = 2 and all iVj = 1, 
restricting the solutions to (at most) quadratic in S and T. 
Example: For (J42j) . we set all Mt = 2, Ni — 1, and continue with 

Ui(S) = a l0 + a a S + a l2 S 2 + VT^S 2 ~(b t0 + b ll S), % = 1,2,3. (44) 

Step ST3: Derive the algebraic system for the coefficients and b^ 

Substituting (|4*UJ) into (j2EJ) gives P(S) + yl — S 2 Q(S), which must vanish iden- 
tically. Hence, equate to zero the coefficients of the power terms in S so that 
P = and Q = 0. 

Example: After substitution of (jBj) into (|42|). the resulting nonlinear algebraic 
system for the coefficients a^- and 6y contains 25 equations (not shown). 
Step ST4: Solve the nonlinear parameterized algebraic system 
In contrast to step S4 we no longer assume that g^m, and 6^ are nonzero (at 
the cost of generating some constant solutions, which we discard later). 
Example: For (|41|). there are 11 solutions. Three are trivial, leading to constant 
Ui. Eight are nontrivial solutions giving the results below. 
Step ST5: Build and test the solitary wave solutions 
Proceed as in step S5. 

Example: The solitary wave solutions of (j41j) are: 

u(x,t) = ±c 2 tanh£, 

v(x, t) = t|c 2 (ci - c 2 ) sech 2 £, (45) 
w(x,t) = -\c 2 {ci - c 2 ) sech 2 £, 

which could have been obtained with the tanh- method of Section [^J 
u(x, t) = ±ic 2 sech£, 

v(x, t) = ±|ic 2 (ci — c 2 ) tanh ^ sech ^, (46) 
w(x,t) = |c 2 (ci - c 2 ) (l - 2sech 2 ^) , 

reported in Gao and Tian (2001); and the two complex solutions 

u(x, t) = ±^ic 2 (sech £ ± i tanh £ ) , 

v(x, t) = |c 2 (ci — c 2 ) sech£ (sech £ ± i tanh£) , (47) 
w(x, t) = — \c 2 (ci — c 2 ) sech£ (sech£ ± i tanh£) . 

In all solutions £ = C\X + c 2 t + A, with Ci,c 2 and A arbitrary. The complex 
conjugates of (}4"Tj) are also solutions. 
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5. Algorithms to compute sn and cn solutions for PDEs 

In this section we give the main steps (labelled CN1-CN5) of our algorithm to 
compute polynomial solutions of (jSJ) in terms of Jacobi's elliptic cosine function 
(cn). Modifications needed for solutions involving the sn function are given at the 
end of this section. Details for steps CN2, CN4 and CN5 are shown in Section |3 
Example: Consider the KdV equation (Ablowitz and Clarkson, 1991), 

u t + auu x + u xxx = 0, (48) 

with real constant a. The KdV equation models, among other things, waves in 

shallow water and ion-acoustic waves in plasmas. 

Step CN1: Transform the PDE into a nonlinear ODE 

Similar to the strategy in Tl and SI, using (Lawden, 1989) 

sn 2 (£; m) = 1 — cn 2 (£; m), dn 2 (^;m) = l-m + mcn 2 ((;m), (49) 

and 

cn'(£; m) = — sn(£; m)dn(£; m), (50) 

one has CN' = -^(1 - CN 2 )(1 - m + m CN 2 ) where CN = cn(£; m) is the 
Jacobi elliptic cosine with argument £ and modulus < m < 1. 
Repeatedly applying the chain rule 

dm dm dCN <9£ f~ ™ t2 x /1 ™ t2 x dm 

W, = dCNlT 4 = ~ c > V d - CN X 1 - m + mCN "Jsraf ■ < 5l > 

system (j2J) is transformed into a nonlinear ODE system. In addition to the Cj, 
the algorithm introduces m as an extra parameter. 
Example: Using (|5ljl to transform (|48j) we have 

(c?(l - 2m + 6m CN 2 ) - c 2 - ac x Ui) U[ , , 

+ 3c?CN(l - 2m + 2mCN 2 )(/; -c?(l - CN 2 )(l - m + mCN 2 )f/f = ^ ' 

Step CN2: Determine the degree of the polynomial solutions 

Follow the strategy in step T2. 

Example: For (gSJ), substituting E/i(CN) = CN Ml into (J52D and equating the 
highest exponents gives l + Mi = — l + 2 M\. Then, M\ = 2, and 

£/i(CN) = a 10 + a n CN + a 12 CN 2 . (53) 

Step CN3: Derive the algebraic system for the coefficients ay- 
Proceed as in step T3. 

Example: For ()48j). after substituting (jS^jl into (|52j). one finds 

an Ci (a a 12 — 2 m c 2 ) = 0, 
ai2 ci (a ai2 — 12 m c 2 ) = 0, 
a n (a a 10 c x — c\ + 2 m c\ + c 2 ) = 0, 
a a 2 x Ci + ai2 (2 a aio ci — 16m c\ — 8 cf + 2 c 2 ) = 0. 
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Step CN4: Solve the nonlinear parameterized algebraic system 

Solve the system as in step T4. 

Example: For c±, c 2 , m, a and a% 2 nonzero, the solution of (|54|) is 



Step CN5: Build and test the solitary wave solutions 

Substitute the results of step CN4 into Reverse step CN1. Test the solutions. 
Example: The cnoidal wave solution of (|4*Hjl is 

u(x,t) = [4ci(l -2m) -c 2 ]/(oci) + (12mc 2 )/(o) cn 2 (ci£ + c 2 t + A; m), (56) 

where c\, c 2 , a, A and modulus m are arbitrary. If any of the parameters in (JHJ) 
are zero, steps CN1-CN5 should be repeated. 
Computation of solutions involving Sn 

To find solutions in terms of Jacobi's sn function, one uses the identities, 

cn 2 (£;m) = 1 - sn 2 (£;m), dn 2 (£; m) = 1 - m sn 2 (£; m), 
sn'(£; m) = cn(£; m) dn(£; m). 

Then, SN' = a/(1 — SN 2 )(1 — mSN 2 ), where SN = sn(^; m) is the Jacobi elliptic 
sine with argument £ and modulus < m < 1. The steps are identical to the cn 
case, except one uses the chain rule 



Since (joTj) and (p)%j) involve roots, as in Sections El and 0] there is no reason to 
restrict the solutions to polynomials in only cn or sn . Solutions involving both 
sn and cn (or combinations with dn) are beyond the scope of this paper. 

Finally, from the sn and cn solutions, sin, cos, sech, and tanh-solutions can be 
obtained by taking the appropriate limits for the modulus (m — > 0, and m — > 1). 
Indeed, sn(£; 0) = sin(£), sn(£; 1) = tanh(f), cn(£; 0) = cos(£), cn(£; 1) = sech(£). 
No need to compute solutions in dn explicitly since cn(-\/m£; 1/m) = dn(£;ra). 

6. Key Algorithms 

In this section we present in a uniform manner the details of steps two, four and 
five of the algorithms in Sections EKD 

6.1. Algorithm to compute the degree of the polynomials 
Step Ml: Substitute the leading-order ansatz 

A tracking variable is attached to each term in the original system of PDEs. Let 
Tr[i] denote the tracking variable of the zth term in 



a w = [Ac\ (1 - 2m) - c 2 ]/(ac 1 ), a n = 0, an = (I2mc\)/a. 



(55) 




(58) 
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F 


R(F) 


T 





S 


1-S' 2 


CN 


(l-CN 2 )(l-m + mCN' A ) 


SN 


(1 - SN*)(1 - mSN') 



Table 1: Values for R(F) in 

The first step of the main algorithms leads to a system of parameterized ODEs 
in U, U', U", . . . , U (m) . These ODEs match the form 

T(F, U(F), U'(F), . . .) + ^/rJf) n(F, U(F), U'(F), . . .) = 0, (59) 
where F is either T, S, CN, or SN, and R(F) is defined in Tabled 
Since the highest degree term only depends on F Mi , it suffices to substitute 

Ui(F) -> F Ml (60) 

into dSH)- 

Example: We use the coupled KdV equations (f2*2*|) as our leading example: 

Tr [1] u t — 6aTr [2] uu x + 2/5Tr [3] vv x — ctTr [4] = 0, 

Tr [5] v t + 3Tr [6] uv x + Tr [7] v xxx = 0. ^ 

Step SI resulted in with n = 0. Substituting (jnOD into (jHU), we get 

(Tr [1] c 2 M 1 - aTr [4] cfM?)^ 1-1 + aTr [4] c?Mi(Mi + l)(Mi + 2)5 Ml+1 

-6aTr [2] c 1 Mi5 ,Mfl - 1 + 2/3Tr [3] c x M 2 S 2M ^ x = 0, 
(Tr [5] c 2 M 2 + Tr [7] c\Ml)S M2 - 1 - Tr [7] c\M 2 (M 2 + 1)(M 2 + 2)S R ' h+1 ' °" ' 

+3Tr[6]c 1 M 2 S Ml+M2 - 1 = 0. 



Step M2: Collect exponents and prune sub-dominant branches 

The balance of highest exponents must come from different terms in (J3J). For each 
equation Aj and for each tracking variable, collect the exponents of F, remove 
duplicates, and non-maximal exponents. For example, Mi — 1 can be removed 
from {Mi + 1, Mi - 1} because Mi + 1 > M x - 1. 

Example: Collecting the exponents of S in (|62|) . we get the unpruned list: 





Ai 




A 2 


Tr[l] : 


{Mi-1} 


Tr[5] : 


{M 2 -l} 


Tr[2] : 


{2Mi-l} 


Tr [6] : 


{Mi + M 2 -1} 


Tr[3] : 


{2M 2 -1} 


Tr [7] : 


{M 2 + 1,M 2 + 1,M 2 + 1,M 2 -1} 


Tr[4] : 


{Mi + 1, Mi + 1, Mi + 1, Mi - 1} 







(63) 
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We prune by removing duplicates and non-maximal expressions, and get 

from Ax : {M x + 1, 1M X - 1, 2M 2 - 1}, 

fromA 2 : {M 2 + 1, Mi + M 2 - 1}. 1 j 



Step M3: Combine expressions and compute relations for Mj 

For each A, separately, equate all possible combinations of two elements. Con- 
struct relations between the Mj by solving for one Mj. 

Example: Combining the expressions in (J64)) . we get 





Ai 




A 2 


Mi 


+ 1 = 2Mi 


- 1 


M 2 + 1 = Mi + M 2 - 1 


Mi 


+ 1 = 2M 2 


- 1 




2Mi 


- 1 = 2M 2 


- 1 





(65) 



We construct relations between the Mj by solving for Mi (in this case): 





Ai 


A 2 


Mi 


= 2 


Mi = 2 


Mi 


= 2M 2 - 2 




Mi 


= M 2 





(66) 



Step M4: Combine relations and solve for exponents Mj 

By combining the lists of expressions in an outer product like fashion, we generate 
all the possible linear equations for Mj. Solving this linear system, we form a list 
of all the possible solutions for Mj. 

Example: Combining the equations in Ai and A 2 , we obtain 

{Mi = 2, Mi = 2}, {Mi = 2, Mi = M 2 }, {M x = 2,M 1 = 2M 2 - 2}. (67) 
Solving, we find 

r Mi = 2 r Mi = 2 

\ M 2 = 2 \ M 2 = Free 1 ' 

Step M5: Discard invalid exponents Mj 

The solutions are substituted into the unpruned list of exponents (in step M2). 
For every solution (without free exponents) we test whether or not there is a 
highest-power balance between at least two different tracking variables. If not, 
the solution is rejected. Non-positive, fractional, and complex exponents are 
discarded (after showing them to the user). Negative exponents (Mj = — pj) and 
fractional exponents (Mj = Pi/qi) indicate that a change of dependent variables 
(u\ = uf Vl or Mj = Mj 1 ^ 1 ) should be attempted in (J3J). Presently, such nonlinear 
transformations are only carried out automatically for single equations. 

Example: Removing the case {Mi = 2, M 2 = Free} from (J68|) . we substitute 
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{Mi = 2, M 2 = 2} into (JHEJ). Leading exponent (3 in this case) occurs for Tr [2] , 
Tr [3] and Tr [4] in Ai, and for Tr [6] and Tr [7] in A 2 . The solution is accepted. 
Step M6: Fix undetermined Mj and generate additional solutions 

When some solutions involve one or more arbitrary Mj we produce candidate 
solutions with a countdown procedure and later reject invalid candidates. 

Based on the outcome of step M5, scan for freedom in one or more of Mj by 
gathering the highest-exponent expressions from the unpruned list in step M2. 
If the dominant expressions are free of any of the Mj, a countdown mechanism 
generates valid integer values for those Mj. These values of Mj must not exceed 
those computed in step M5. Candidate solutions are tested (and rejected, if 
necessary) by the procedure in step M5. 

Example: The dominant expressions from IjfiUjl with {Mi = 2, M 2 = 2} are 





Ai 






A 2 


Tr[2] 


{2Mi 


-1} 


Tr[6] : 


{Mi + M 2 - 1} 


Tr[3] 


{2M 2 


-1} 


Tr[7] : 


{M 2 + 1} 


Tr[4] 


{Mi 4 


-1} 







(69) 



Substituting Mi = 2, the highest exponent (3 in this case) matches for Tr[2] 
and Tr[4] in Ai when M 2 < 2. The highest exponent (M 2 + 1) matches for 
Tr[6] and Tr[7] in A 2 . 

A countdown mechanism then generates the following list of candidates 



M 1 = 1 f Mi = 1 j Mi = 2 f Mi = 2 

M 2 = 1 \ M 2 = 2 \ M 2 = 1 \ M 2 = 2 

Verifying these candidate solutions, we are left with 

Mi = 2 f M 2 

M 2 = 1 1 M 2 = 2 



(70) 



(71) 



Notice that for the new solution {Mi = 2, M 2 = 1} only the exponents corre- 
sponding to Tr [2] and Tr [4] in Ai are equal. 

Currently, for the mixed tanh-sech method, the code sets Mj = 2 and iVj = 1. 



6.2. Algorithm to analyze and solve nonlinear algebraic systems 

In this section, we detail our algorithm to analyze and solve nonlinear parameter- 
ized algebraic systems (as generated in step 3 of the main algorithms). Our solver 
is custom-designed for systems that are (initially) polynomial in the primary un- 
knowns (cijj), the secondary unknowns (cj), and parameters (m, a, /3, 7, . . .). 

The goal is to compute the coefficients in terms of the wave numbers 
Cj and the parameters m, a, /3, etc. In turn, the Cj must be solved in terms of 
these parameters. Possible compatibility conditions for the parameters (relations 
amongst them or specific values for them) must be added to the solutions. 

Algebraic systems are solved recursively, starting with the simplest equation, 
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and continually back-substituting solutions. This process is repeated until the 
system is completely solved. 

To guide the recursive process, we designed functions to (i) factor, split, and 
simplify the equations; (ii) sort the equations according to their complexity; 
(iii) solve the equations for sorted unknowns; (iv) substitute solutions into the 
remaining equations; and (v) collect the solution branches and constraints. 

This strategy is similar to what one would do by hand. If there are numerous 
parameters in the system or if it is of high degree, there is no guarantee that our 
solver will return a suitable result, let alone a complete result. 
Step Rl: Split and simplify each equation 

For all but the mixed tanh-sech algorithm, we assume that the coefficients a^Mi 
of the highest power terms are nonzero and that q, m, a, /3, etc. are nonzero. For 
the mixed sech-tanh method, a,^ = a>i2 and 6^ = hi are allowed to be zero. 

We first factor equations and set admissible factors equal to zero (after clearing 
possible exponents). For example, {0i</>2 3 03 2 = 0} — » {fa = 0, fa = 0,03 = 0}, 
where fa is a polynomial in primary and secondary unknowns along with the 
parameters. Equations where non-zero expressions are set to zero are disgarded. 

Example: Consider (J34)) . which was derived in the search for sech-solutions of 
for case Mi = M 2 = 2. Taking di 2 , a 22 , Cj, c 2 , a, /?, to be nonzero, splitting 
equations, and removing non-zero factors leads to 

a 12 - Ac\ 

a 2 i = V (3a 10 ci + c\ + c 2 ) 
ai2G2i + 2ana 22 — 2a 2 ic\ 
3aanai2 — j3a 2 ia 22 — aauc\ 
?>aa\ 2 — /3a 22 — Qaa\ 2 c\ 
6aaioanCi — 2/3a 2 oa 2 iCi + aauc\ — anc 2 
3a n a 21 ci + 6ai a 22 ci + Sa 22 c\ + 2a 2 2C2 
^aa^ci + 6aaioai2Ci — /3a 21 ci — 2/5a 2 oa22Ci + Aaa^cf — a\ 2 c 2 

where V is the logical or. 

Step R2: Sort equations according to complexity 

A heuristic measure of complexity is assigned to each fa by computing a weighted 
sum of the degrees of nonlinearity in the primary and secondary unknowns, pa- 
rameters, and the length of fa. Linear and quasi-linear equations (with products 
like axia 2 i) are of lower complexity than polynomial equations of higher degree 
or non-polynomial equations. Solving the equation of the lowest complexity first, 
forestalls branching, avoids expression swell, and conserves memory. 



0. 
0. 
0. 
0. 

o, 
o, 
o, 

0. 



(72) 
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Example: Sorting (|72jl . we get 

a 12 - 4df 

3aanai2 — Pa 2 ±a 22 — aancf 
«i2G2i + 2ana 2 2 — 2a 2 ic\ 
a 21 = V (3ai ci + + c 2 ) 
3aa^ 2 — /?«22 — 6aai2C^ 
6aaioanCi — 2(3a 2 oa 2 iC\ + aanCj — a\\c 2 
3aua 21 Ci + Qa w a 22 Ci + %a 22 c\ + 2a 2 2C2 
3aa^ci + 6aaioai2Ci — /3a2i c i — 2/3a2ofl22Ci + 4:aai 2 cf — a\ 2 c 2 

Step R3: Solve equations for ordered unknowns 

The ordering of unknowns is of paramount importance. The unknowns from 
the first equation from step R2 are ordered so that the lowest exponent primary- 
unknowns precede the primary- unknowns that the equation is not polynomial in. 
If there are not any primary-unknowns, the lowest exponent secondary-unknowns 
precede the secondary-unknowns that the equation is not polynomial in. Like- 
wise, in the absence of primary- or secondary- unknowns, the lowest exponent 
parameters precede the non-polynomial parameters. 

The equation is solved using the built-in Mathematica function Reduce, which 
produces a list of solutions and constraints. Constraints of the form a ^ b (where 
neither a or b is zero) are pruned, and the remaining constraints and solutions 
are collected. 

Example: In this example, a\ 2 — Ac\ = is solved for a\ 2 and the solution 
a 12 = Ac\ is added to a list of solutions. 
Step R4: Recursively solve the entire system 

The solutions and constraints from step R3 are applied and added to the previ- 
ously found solutions and constraints. In turn, all the solutions are then applied 
to the remaining equations. The updated system is simplified by clearing com- 
mon denominators in each equation and continuing with the numerators. Steps 
R1-R4 are then repeated on the simplified system. 

Example: Substituting a\ 2 = 4c^ and clearing denominators, we obtain 

Pa 2 \a 22 — llaaiicl = 

a ll a 22 + a 21 c l = 0, 

a 21 = V (3ai ci + c\ + c 2 ) = 0, 

f3a 2 22 - 24a4 = 0, (74) 

6aaioanCi — 2(3a 2 Qa 2 \Ci + aancf — a\\c 2 = 0, 

3ana 2 iCi + 6a 10 a 22 ci + %a 22 c\ + 2a 2 2C2 = 0, 

3aa^ — (3a\ x — 2Pa 20 a 22 + 24aa w cl + lQac\ — 4ciC 2 = 0. 

The recursive process terminates when the system is completely solved. The 
solutions (including possible constraints) are returned. 



0. 

o, 
o, 
o, 
o, 
o. 
o, 

0. 



(73) 
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Repeating steps R1-R4 seven more times the global solution of (pMjl is obtained: 
a w = -(4c? + c 2 )/(3ci), a u = 0, a X2 = 4:cf, 

(75) 

a 20 = ±(4ac? + (1 + 2a)c 2 )/(c 1 y/6a(3), a 2 i = 0, a 22 = =F2c^ ^/6a//?. 

where ci, c 2 , a and /3 are arbitrary. 

This solution of corresponds to the Mi = 2, M 2 = 1 case given in (j3SJ). 

6.3. Algorithm to build and test solutions 

The solutions to the algebraic system found in Section 16.21 are substituted into 

Mi Ni 

u,(x) = J>^'(£) + V^Wj^i^O. (76) 

j=0 j=0 

where F and -R(-F) are defined in Section |6~T1 The constraints on the parameters 
(m,a,/3, etc.) are also collected and applied to system (J3j). 

Since the algorithm used to solve the nonlinear algebraic system continually 
clears denominators, it is important to test the final solutions for Uj. While 
Mathematica's Reduce function generates constraints that should prevent any 
undetermined or infinite coefficients a^- after back-substitution, it is still prudent 
to check the solutions. 

To present solutions in the simplest format, we assume that all parameters 
(ci,m,a, j3, etc.) are positive, real numbers. This allows us to repeatedly apply 
rules like \fa^ — > a, y/—a 2 — > ia, \f—fi — > i \ffi and \f—{c\ + c 2 ) 2 — > i (ci + c 2 ). 

We allow for two flavors of testing: a numeric test for complicated solutions and 
a symbolic test which guarantees the solution. In either test, we substitute the 
solutions into (jSJ) after casting the solutions into exponential form, i.e., tanh£ — ► 
(e* - e - *)/(e* + e"«) and sech^ -> 2/(e« + e"«). 
For the numeric test of solutions: 

• after substituting the solution, substitute random real numbers in [—1, 1] 
for Xi,Ci, and A in the left hand side of © , 

• expand and factor the remaining expressions, 

• substitute random real numbers in [—1, 1] for arbitrary a^-, bij,m,a, . . . , 

• expand and factor the remaining expressions, 

• if the absolute value of each of the expressions < e m MachinePrecision/2, 
then accept the solution as valid, else reject the solution (after showing it 
to the user). 

Mathematica evaluates ya? — > a when a is numeric, but does not evaluate vx — > 
\a\ when a is symbolic. Our simplification routines use \/~a? = a instead of \a\ 
when a is symbolic. This has two consequences: (i) valid solutions may be missed, 
and (ii) solutions have a 1/2 probability of evaluating to matching signs during 
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the numeric test. The numeric test being inconclusive, we perform a symbolic 
test. 

For the symbolic test of solutions: 

• after substituting the solution, expand and factor the left hand side of (jHJ), 

• apply simplification rules like \fa? —>■ a, \J —a 2 — > i a, 1 — sech 2 £ — ► tanh 2 £, 
and sn 2 (x; m) — > 1 — cn 2 (x; m), 

• repeat the above simplifications until the expressions are static, 

• if the final expressions are identically equal to zero, then accept the solution, 
else reject the solution and report the unresolved expressions to the user. 

7. Examples of solitary wave solutions for ODEs and PDEs 

The algorithms from SectionsEKJwere implemented in our Mathematica package 
PDESpecialSolutions .m, which was used to solve the equations in this Section. 

7.1. The Zakharov-Kuznetsov KdV-type equations 

The KdV-Zakharov-Kuznetsov (KdV-ZK) equation, 

u t + auu x + u xxx + u xyy + u xzz = 0, (77) 

models ion-acoustic waves in magnetized mult i- component plasmas including 
negative ions (see e.g. Das and Verheest, 1989). 

With PDESpecialSolutions (Tanh and Sech options) we found the solution 



8ci(c? + c\ + cl) - c 4 12(c? + c 2 . + c 2 ;, 

u(x, y, z, t) = 2 ^ ^ 2 - ^ tanh 2 £, 

aci a 



4ci(c? + c 2 2 + cl) + c 4 + 12(c 2 + c 2 2 + c 2 ) gech ^ 



(78) 



olc\ a 



where £ = c\x + c 2 y + c 3 z + c^t + A, with c\, c 2 , c 3 , c 4 , A and a arbitrary. 
For c 2 = c 3 = and replacing C4 by c 2 , one gets the solitary wave solution 



. . 8ci — c 2 12ci 2/ \ 

t) = tanh (cix + c 2 t + A), 

aci a 

4cl + c 2 12c 2 



, , 12c 2 (79) 
H sech 2 (cix + c 2 t + A), 



aci a 



of the ubiquitous KdV equation 

The function PDESpecialSolutions does not take boundary or initial condi- 
tions as input. One can a posteriori impose conditions on solutions. For example, 
requiring linx^-too u(x, t) = in (fTH|) would fix c 2 = — 4cf. 

For the modified KdV-ZK equation (Das and Verheest, 1989), 

u t + au 2 u x + u xxx + u xyy + m K2 = 0, (80) 
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using the Tanh and Sech options, PDESpecialSolutions returns 



u(x,y, z,t) = ±i y 6(cf + c 2 + c§)/a!tanh£, (81) 
with £ = C\X + c 2 y + c 3 z + 2c\{c\ + c\ + c£)t + A, and 



u(x, y, z, t) = ±y 6(cf + c| + c|)/a sech £, (82) 

with £ = ax + c 2 y + c 3 z - ci(c{ + c 2 2 + dj)t + A. For c 2 = c 3 = 0, flS3} and 
reduce to the well-known solitary wave solutions 



u(x,t) = ±ic 1 ^/6/atanh(cix + 2c\t + A), (83) 
u(x,t) = ±ciA/6/asech(ciX - c\t + A) (84) 

(ci,A and a arbitrary real numbers) of the modified KdV (mKdV) equation 
(Ablowitz and Clarkson, 1991), 

u t + au 2 u x + u xxx = 0. (85) 
For a three dimensional modified KdV (3D-mKdV) equation, 

u t + 6u 2 u x + u xyz = 0, (86) 
one obtains the solitary wave solution 

u(x, y, z, t) = ±y/c 2 c 3 sech(ci x + c 2 y + c 3 z - c x c 2 c z t + A), (87) 
where c\, c 2 , c 3 and A are arbitrary. 

7.2. The generalized Kuramoto-Sivashinsky equation 

Consider the generalized Kuramoto-Sivashinsky (KS) equation (see e.g. Parkes 
and Duffy, 1996): 

Ut + uu x + u xx + au xxx + u xxxx = 0. (88) 

Ignoring complex solutions, PDESpecialSolutions (Tanh option) automatically 
determines the special values of the real parameter a and the corresponding 
closed form solutions. For a = 4, 

u(x, t) = 9 ± 2c 2 ± 15 tanh£ - 15 tanh 2 £ =F 15tanh 3 £, (89) 

with £ = =f|x + c 2 t + A. For a = 

, . 45^4418c 2 45 , ^ 45 , 2> 15 , ^ , . 

u(x, t) = — — ± — tanh £ — tanhV ± = tanh 3 £, (90) 

47^ 47^47 47^ 47^ 
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where £ = ±^jp + c 2 t + A. For a = 



, , 2(30 =F 5329ca) 75 , ^ 60 , a> 15 , , , 

= — — = ± = tanh£ — tanh¥± — tanhY, (91) 

73^ 73^ 73^ 73^ 



where £ = ±^g£ + c 2 t + A. 

The remaining solutions produced by PDESpecialSolutions are either com- 
plex (not shown here) or can be obtained from the solutions above via the in- 
version symmetry of IjBBJl : u — > —u, x — > —x, a — > —a. 

A separate run of the code after setting a = in (jHHj) yields 



/ \ „ /19 135 /ll , „ 165 /ll , 

u(x,t) = -2\ — c 2 \ — tanh£ + \ — tanh 3 £, (92) 

v ' ; V 11 19 V 19 S 19 V 19 S ' V ' 

with £ = \\J^x + c 2 t + A. In all solutions above c 2 is arbitrary. 
7.3. Coupled KdV equations 

In Sectional we gave the sech-solutions for the Hirota-Satsuma system (|22[) . Here 
we list the tanh, cn and sn solutions for (|22|) computed by PDESpecialSolutions 
(Tanh, JacobiCN and JacobiSN options): 

u(x, t) = 2Cl - ° 2 - 2c? tanh 2 (£), 
3ci 

v(x, t) = ±J [8acf + 2(1 + 2a)ciC 2 ]//3 tanh(£), 

8c 3 - co (93) 

M (*,t) = ^-^-4c 2 tanh 2 (0, 
3ci 

*(*,*) = ± 8 ^ ~ ( L±f )C2 T 2c 2 v^tanh 2 (0, 



u(x,t) = (1+ ^ )C ' C2 -2mc 2 sn 2 (e;m), 
3ci 

w(x, t) = ±a/ [4am(l + m)c\ + 2(1 + 2a)m Cic 2 }/ (3 sn(f , m), 

4(1 + m)c\ - c 2 2 2 (94) 

r) = 4m c : sn (f ; m), 

3ci 

4c* (1 + m)c? — (1 + 2a)c 2 9 , —- 9 .. 

u(z, t = ±— ^ ^=== ^ T 2c 2 v^oT^ sn 2 (£; m), 

ci v6ap 
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= (1 2m )4 02 +2 mclcn 2 (Z;m), 
3ci 

t) = ±\J[4am(2m - l)c\ - 2(1 + 2a)mciC 2 ]/7? cn(£; m), 

, , 4(l-2m)c?-c 2 „ o 2 / a \ ^ 
3ci 

4a(l — 2m)c? — (1 + 2a)c 2 9 , -r-z 0/J . 

v(x, t) = ±—± ^ ± 2c? v^T^ cn 2 (£; m , 

with £ = cix + c 2 t + A, and ci, c 2 , a, /3, A, and modulus m arbitrary. These 
solutions correspond with those given in Fan and Hon (2002). 

With the SechTanh option we obtained two dozen (real and complex) solu- 
tions. The real solutions coincide with the ones given above. 

Another coupled system of KdV-type equations was studied by Guha-Roy 
(1987) 

Ut + avv x + (3uu x + ^u xxx = 0, 

(96) 

v t + 5(uv) x + evv x = 0, 

where a through e are real constants. The package PDESpecialSolutions (Sech 
option) computed: 

, , 4e 2 7c? + (4aS + e 2 )c 2 12e 2 7c? , ,, 

u(x, t) = \ '— + sech 2 (c x x + c 2 t + A , 

Act A ' 

v(x t) = wuri+y-mcA _ « secm (9 f 

where £ = Cix + c 2 t + A, A = 4a<5 2 + /?e 2 , with ci, c 2 , A and a through e arbitrary. 
For e = 0, (J9H)) reduces to Kawamoto's system; for e = 0, 5 = —2 to Ito's system. 
Neither of these systems has polynomial solutions in sech or tanh . 
7.4. The Fisher and FitzHugh-Nagumo equations 

For the Fisher equation (Malfhet, 1992), 

u t - u xx - u(l - u) = 0, (98) 

with PDESpecialSolutions (Tanh option) we found the (real) solution 

u(x,t) = \ ± |tanh£ + ^tanh 2 £, (99) 

with £ = ±tA?^ ± j^t + A. In addition, there are 4 complex solutions. 

Obviously, PDESpecialSolutions handles ODEs also. For example, we can 
put the FitzHugh-Nagumo (FHN) equation (Hereman, 1990), 

Ut — u xx + u(l — u)(a — u) = 0, (100) 
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where — 1 < a < 1, into a travelling frame, 

pv s + V2v zz -V2v(l-V2v)(a-V2v) = 0, (101) 

with v(z) = v(x — ^=t) = \/2u(x,t). Ignoring the inversion symmetry z — > 
—z, j3 — > —f3 of (jlOljl . we find with PDESpecialSolutions (Tanh option) 



v(z) 



2V2 



(3 + (/3-2) tanh[^-(2 - p)z + A] 



(102) 



if a — (3 — 1; 



v(z) 



(£+3 
2V2 



\[2 

1 - tanh[-^— (/3 + 2)z + A] 



if a = P + 2; and 



v/2 

viz) = — = 1 + tanhf— z + A] 
2y/2 1 4 



(103) 



(104) 



if a = =(/3+l). In these solutions (see e.g. Hereman, 1990) P and A are arbitrary. 



7.5. A degenerate Hamiltonian system 

Gao and Tian (2001) considered the following degenerate Hamiltonian system, 



u t — u x — 2v = 0, 

v t — 2euv = 0, e = ±1, 



(105) 



which was shown to be completely integrable by admitting infinitely many con- 
served densities. Our code does not find sech-solutions. With the SechTanh op- 
tion, PDESpecialSolutions returns the solutions: 



u(x, t) = — ec 2 tanh£, 

v(x,t) = |ec 2 (ci -c 2 )sech 2 £, 

which could have been obtained with the tanh-method in Section and 

u(x, t) = \1c2e (sech £ + i tanh £), 

v(x,t) = |c 2 (ci — c 2 )e sech£ (sech£ + ztanh£), 



(106) 



(107) 



plus their two complex conjugates. There are no constraints on ci, c 2 , and e, and 
£ = C\X + c 2 t + A. The above solutions were reported in Gao and Tian (2001). 
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7.6. The combined KdV-mKdV equation 

The combined KdV-mKdV equation (see Gao and Tian, 2001) 

u t + 6auu x + 6j3u 2 u x + r yu xxx = 0, (108) 

describes a variety of wave phenomena in plasma, solid state, and quantum 
physics. We chose this example to show that ODEs of type (|27|). which are free 
of \/l — S 2 , can admit mixed tanh-sech solutions. 

First, PDESpecialSolutions with the Tanh option, produces 



u(x,t) = -— ±y^citanh(cix + ^(3^ + 4/3 7 cf)t + A). (109) 

Next, with the Sech option, PDESpecialSolutions computes 

u(x, t) = -^± y| Cl sech[ Cl x + ^(3a 2 - 2{3 1C \)t + A]. (110) 

Third, with the SechTanh option, PDESpecialSolutions finds 

ot i — 

u(x,t) = -— + |y^ci (sech^ iitanhO, (HI) 

and 

Qt I — 

u(x,t) = -— - ^y^ci (sech ^=Fi tanh 0, (112) 

where ^ = C\X + |^-(3a 2 + /?7cf)t + A. In all solutions ci, A,a,/3 and 7 are 
arbitrary. The solutions were reported in Gao and Tian (2001), although there 
were minor misprints. 
7.7. The Duffing Equation 

Duffing's equation (Lawden, 1989), 

u" + u + cm 3 = 0, (113) 
models a nonlinear spring problem. Its cn and sn solutions 



2m ex 

uix) = ±a / r — cnf . + A: m), e = ±1, 

W (l-2m)a VI -2m ; 



u(x) = ±\l — 2m sn( eX + A;m), e = ±1, 



(114) 



are computed by PDESpecialSolutions with the JacobiCN and JacobiSN op- 
tions. There are four sign combinations in (|114|) . Since < m < 1, the cn solution 
is real when a > and m < |. The sn solution is real for a < 0. Such conditions 
are not automatically generated. During simplifications the code assumes a > 
(see Section HT21 for details). 

Initial conditions fix the modulus in (|114|) . For example, u(0) = a and u(0) = 
lead to u{x) = a cn(vl + aa 2 x; (aa 2 )/(2 + 2aa 2 )). 
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7.8. A class of fifth-order PDEs with three parameters 

To illustrate the limitations of PDESpecialSolutions consider the family of 
fifth-order KdV equations (G6kta§ and Hereman, 1997), 

u t + au 2 u x + (3u x u xx + "yuu xxx + u xxxxx = 0, (115) 

where a, /3, and 7 are nonzero parameters. 

An investigation of the scaling properties of ()115j) reveals that only the ratios 
and - are important, but let us proceed with (|115|) . 

Special cases 

Several special cases of ()115|) are well known (for references see G6kta§ and 
Hereman, 1997). Indeed, for a = 30, (3 = 20, 7 = 10, Eq. f)115j) reduces to 

u t + 30u 2 u x + 20u x u xx + 10uu xxx + u xxxxx = 0, (H6) 

which belongs to the completely integrable hierarchy of higher-order KdV equa- 
tions constructed by Lax. Equation ()116|) has two tanh-solutions: 

u(x,t) = 4c 2 - 6c 2 tanh 2 (ci£ - h<oc\t + A), (117) 

and 

u(x, t) = aio - 2c\ tanh 2 [cia; - 2(15a 2 ci - A$a w c\ + 2%c\)t + A], (118) 

where a 10 , C\, A are arbitrary. 

For a = (5 = 7 = 5, one gets the equation, 

u t + 5u 2 u x + hu x u xx + huu xxx + u xxxxx = 0, (119) 

due to Sawada and Kotera (SK) and Dodd and Gibbon, which has tanh-solutions: 

u(x, t) = 8cj - 12c 2 tanh 2 (cia; - \Qc\t + A), (120) 

and 

u(x, t) = 010 - 6c 2 tanh 2 [cix - (5a 2 ci - A0a 10 cl + 76cf)t + A], (121) 

where aio, ci, A are arbitrary. 

The KK equation due to Kaup and Kupershmidt, 

Ut + 20u 2 u x + 25u x u xx + 10uu xxx + u xxxxx = 0, (122) 

corresponding to a = 20,(3 = 25,7 = 10, and again admits two tanh-solutions: 

3 

u(x,t) =c\- -c\ tanh 2 (c x x - c\t + A) , (123) 

and 

u(x,t) = 8cj - 12c 2 tanh 2 ( Cl x - 176cft + A), (124) 
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with Ci, A arbitrary, but no additional arbitrary coefficients. 
The equation 

u t + 2u 2 u x + Qu x u xx + 3uu xxx + u xxxxx = 0, (125) 

for a = 2, (3 = 6, 7 = 3, was studied by Ito. It has one tanh solution: 

u(x, t) = 20c? - 30c? tanh 2 (cix - 96cft + A), (126) 

again with c% and A arbitrary. PDESpecialSolutions (Tanh option) produces 
all these solutions. 

General case 

Eq. ()115|) is hard to analyze by hand or with the computer. After a considerably 
amount of time, PDESpecialSolutions (Tanh option) produced the solutions 
given below (but not in as nice a form). Our write-up of the solutions is the 
result of additional interactive work with Mathematica. 
The coefficients a 10 , an, and a i2 in 

u(x, t) = aio + a n tanh(£) + ai 2 tanh 2 (£), (127) 

with £ = C1X + C2 + A, must satisfy the following nonlinear algebraic system with 
parameters C\ y C2, a, (3, and 7 : 

aa\ 2 + 6j3ai 2 cl + 127a 12 c? + 360c^ = 0, 
an {aa\ 2 + 2(3a 12 cj + 67012c 2 + 24c^) = 0, 
an (aa 2 ci - 270^ + 2(3a u cl + 16cf + c 2 ) = 0, 
an (aa n + 6aai ai 2 + 67010c 2 — 12/3ai 2 c 2 — 187ai 2 c? — 120c^) = 0, 

2aa 2 1 ai 2 + 2aa w a\ 2 + /?a 2 x c? + 37a 2 1 dj ! + 127Oi ai 2 c 2 (^8) 

-%f3a 2 l2 c\ - 8-fa 2 l2 cl - 480ai 2 c^ = 0, 
aaioa 11 ci + aa 10 ai 2 ci — pa 11 c 1 — 7a 11 c 1 — 87aioai 2 c 1 + 2pa 12 c 1 

+ 136ai 2 c^ + a i2 c 2 = 0. 

Assuming nonzero ai 2 , c\, c 2 , a, /3, and 7, two cases must be distinguished: 
Case 1: an = 0. In turn, this case splits into two sub-cases: 
Case la: 

3 

an = 0, a 12 = --ai , c 2 = c\(2Ac\ - Pa w ), (129) 
where aio must be one of the roots of 

aa 2 l0 - 4f3a l0 cl - 87010c 2 + 160cf = 0. (130) 

Case lb: 

an = 0, a 12 = c?, c 2 = (a 2 a? Ci — 8a7ai c? + 16ac? + 127 2 c?), (131) 

a a 
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provided that 

P = 

7 

Case 2: a n ^ 0. Then 



p = I(i0a- 7 2 ). (132) 



01 = 3^2 (8/?2 + 38/?7 + 3972) ' ai2 = _ 2^^ C2 ' (133) 
provided (3 is one of the roots of 

(104/3 2 + 886/37 + 14877 2 ) (520/5 3 + 2158/3% - 1103/37 2 - 88717 3 ) = 0. (134) 

Thus, case 2 also splits into two sub-cases: 
Case 2a: If (3 2 = -^(886^7 + 1487 7 2 ), then 



1 52(4378/j + 9983 7 2 , 336 2 

a = (2p + 07 7, am = : rCi, an = ±— c l5 

26 V M 77(958/3 + 22137) x ' 2/3 + 37 

168 2 _364(1634/3 + 38517) B 

° 12 ~ ~2/3 + 37 Cl ' ° 2 ~ 2946/3 + 67157 Cl ' 

where (3 is any root of 104/3 2 + 886/37 + 14877 2 = 0. 
Case 2b: If /3 3 = ^(1103/?7 2 + 887I7 3 - 2158/? 2 7), then 

1 <*fP^-i*R 2n 28(1066/3 2 + 5529/37 + 6483 7 2 ) 2 

a = — — -(8/3 + 38/37 + 397 ), aio 



392 v ^ ^' ' " (2/3 + 37)(6/3 + 237) (26/3 + 8I7) 



(135) 



2 _ 28224(26/3- 177) (4/3 -7) 4 168 2 

ttn ~ (2/3 + 3 7 ) 2 (6/3 + 23 7 ) (26/3 + 81 7 ) Cl ' ° 12 ~ ~2/3 + 3 7 Cl ' (136) 

8(188900114/3 2 + 1161063881/37 + 17922619777 2 ) 5 
° 2 ~ 105176786/3 2 + 632954969/37 + 9598334737 2 

where (3 is any root of 520/3 3 + 2158/3 2 7 - 1103/37 2 - 887I7 3 = 0. 
8. Other Algorithms and Related Software 

8.1. Other perspectives and potential generalizations 

The algorithms presented in this article can be extended in several ways. For 
instance, one could modify the chain rule in step Tl (SI, ST1, or CN1) to 
compute other types of solutions or consider more complicated polynomials than 
those used in step T2 (S2, ST2, or CN2). Both options could be used together. 

With respect to the first option, it suffices to know the underlying first-order 
differential equation of the desired fundamental function in the polynomial so- 
lution. Table El summarizes some of the more obvious choices. 



D. Baldwin et al: Symbolic computation of Tanh, Sech, Cn, and Sn Solutions 30 



Function 


Symbol 


ODE (j/' = f) 


Chain Rule 


tanh(£) 


T 


/ i 2 

v =1 - V 


^l =c ,(1_ t2 )^ 


sech(£) 


S 


y^-yV 1 - y 2 


^ = - c ,sVl-s^ 


tan(£) 


TAN 


y' = l + y 2 


^ = Ci (l + TAN 2 )^ 


exp(C) 


E 


y'=y 




cn(£; to) 


CN 


y' = - v /(l-y 2 )(l-m+my 2 ) 


-L = ^. v /(l_ CN 2 )(1 _ TO+mCN 2 ) _g_ 


sn(£; to) 


SN 


y' = y/(\ - y 2 )(l - my 2 ) 


^=C,V(1-SN 2 )(1-TOSN 2 )^ 




P 


y' = ±V^y 3 ~ 92y- 93 


£-=±c^4y3-g 2 y- g3 $ 



Table 2: Functions with corresponding ODEs and chain rules. V(x\ g 2 , gs) is the Weierstrass 
function with invariants g 2 and 173. 



Several researchers, including Fan (2002abc) and Gao and Tian (2001), seek 
solutions of the form 

Mi 

Ui{x,t) = Ui(t)=J2 a iM0 j i £ = c 1 x + c 2 t + A, (137) 
i=i 

where w(£) is constrained by a Riccati equation, 

w'(£) = b + ew 2 (^), e = ±1, b real constant. (138) 

Ignoring rational solutions, (|138|) has the following solutions 

= a tanh(a£ + c), if e = —1, b = a 2 , 

w(£) = a coth(a£ + c), if e = —1, b = a 2 , 

w(£) = a tan(a£ + c), if e = 1, b = a 2 , ^ 

u>(£) = a cot(a£ + c), if e = —1, 6 = —a 2 . 

So, (|137|) is polynomial in tanh£, tan£, coth£, or cot£. The integration constant 
c gets absorbed in A, and the constant a (or b) is an extra parameter in the 
nonlinear algebraic system for the a^. For single PDEs, Yao and Li (2002ab) 
consider solutions of the form 

M M 

u(x, t) = U(0 = *M0 j + E h J^O<r(CV ; - (140) 

j=0 j=0 

where w(£) and z(£) satisfy the Riccati equations 

™'(0 = -™(0*(0, ^(0 = 1-^(0- ( 141 ) 

Since w(£) = sech(£),z(£) = tanh(£) this approach is similar to the sech-tanh 
method in Section 0J 
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Generalizing further, Fan (2002b, 2003abc), Fan and Hon (2002, 2003a), and 
Hon and Fan (2004b) take 

y'(0 = V b o + hy + b 2 y 2 + b 3 y 3 + hy 4 , b { constant, (142) 

which covers the functions sech, sec, tanh, tan, cn, sn, and V. The parameters hi 
are added to the nonlinear algebraic system, which makes such systems hard 
to solve without human intervention. Most often, such complicated nonlinear 
algebraic systems are solved interactively with the aid of Mathematica or Maple. 
To avoid unmanageable systems, Mj(< 2) is often fixed in ()137j) . Chen and 
Zhang (2003ab), Fan and Dai (2003), and Sirendaoreji (2003, 2004) use variants 
of ()142|) to compute polynomial and rational solutions in terms of tanh, sech, tan, 
Jacobi's elliptic functions, etc. 

Zheng et al. (2002) introduce a clever method to compute mixed tanh-sech 
solutions for the combined KdV-Burgers equations. They seek formal solutions, 

M M 

u(x,t) = C/(0 = a + ^2bjSm j w(£) + cosw(£) smw(£) j ~\ (143) 

j=i j=i 

subject to ^| = sinu;(£) which, upon integration, gives sinu;(£) = sech(£) and 
cos«;(£) = ±tanh(£). Alternatively, one can use ^| = cosu>(£), which leads to 
cosw;(£) = — sech(£) and sinu>(£) = ±tanh(£). 
Liu and Li (2002a) seek solutions of the forms 

M M M 

u(o = E a i s <o j , u(o = e *3 H0 j + E h i cn ^) 

j=0 j=0 j=0 

M M M 

U{£) = E % sn (0 j + E A i cn & sn (0 i_1 + E b > dn (^ sn (0 i_1 (144) 

3=0 3=0 3=0 

M 

+E^ cn (o dn (o H0 j ~ 2 , 

3=0 

which generalize the Jacobi elliptic function method in Section 0] 
With respect to the second option, Gao and Tian (2001) consider 

Mi Ni 

Ui(x,t) = aij(x, t) tanh%(x, t) + E K( x i f ) sech^(x, t) tanh j ^(x, t), (145) 

3=0 3=0 

where ^f(x, t) is not necessarily linear in x and/or t. Of course, ()145j) arises from 
recasting the terms in ()39)1 in a slightly different way than (|40|). Restricted to 
travelling waves, \l/(x, t) = C\X + c^t + A, both forms are equivalent. 

Our algorithms could be generalized in many ways. With considerable effort, 
solutions involving complex exponentials multiplied by tanh or sech functions 
could be attempted. A solution to the nonlinear Schrodinger equation is of this 
form. Fan and Hon (2003b), Hon and Fan (2004a) and Fan (2003bc) give exam- 
ples of complex as well as transcendental equations solved with the tanh-method. 
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8.2. Review of Symbolic Algorithms and Software 

There is a variety of methods to find solitary wave solutions and soliton solu- 
tions of special nonlinear PDEs. See e.g. Hereman and Takaoka (1990), Estevez 
and Gordoa (1995, 1998), and Helal (2002). Some of these methods are straight- 
forward to implement in computer algebra systems (CAS). 

The most comprehensive methods of finding exact solutions for ODEs and 
PDEs are based on similarity reductions via Lie point symmetry methods. These 
methods are hard to fully automate (for publications and software see e.g. 
Cantwell, 2002; Hereman, 1996; Hydon, 2000). Most CAS have tools to solve 
a subset of linear and nonlinear PDEs. For example, Mathematica's DSolve can 
find general solutions for linear and weakly nonlinear PDEs. Available within 
MuPAD, the code pdesolve uses the method of characteristics to solve quasi- 
linear first order PDEs. Maple offers the packages ODEtools (for solving ODEs 
using classification, integrating factor and symmetry methods) and PDEtools, 
which contains the function pdesolve to find exact solutions of some classes of 
PDEs. For information consult Cheb-Terrab (1995, 2001ab). 

The methods presented in this paper are different from these efforts. Our 
algorithms and software only compute specific solutions of nonlinear PDEs which 
model travelling waves in terms of the tanh, sech, sn and cn functions. Our code 
can handle systems of ODEs and PDEs with undetermined parameters. 

To our knowledge, only four software packages are similar to ours. The first 
package is ATFM by Parkes and Duffy (1996), who automated to some degree the 
tanh-method using Mathematica. In contrast to ATFM, our software performs 
the computations automatically from start to finish without human intervention. 
In our code, the number of independent variables X{ is not limited to one space 
variable x and time t; our code handles any number of dependent variables. 

The second package is RATH by Li and Liu (2002), which automates the tanh- 
method. In contrast to our code, RATH only works for single PDEs. Extensions 
to cover systems of PDEs and sech solutions are under development. Surpassing 
our code, RATH can solve PDEs with an unspecified degree of nonlinearity and 
deal with negative and fractional exponents. 

Table El compares the performance of PDESpecialSolutions .m and RATH. 
The solution times are comparable, yet occasionally there is a mismatch in the 
number of solutions computed. This is due in part to the representation of solu- 
tions. Occasionally special solutions are generated although-after inspection by 
hand-they are included in more general solutions. 

Liu and Li (2002a) present the Maple code AJFM to automate the Jacobi elliptic 
function method for single PDEs. This package seeks solutions of the form (|144|) . 

The codes RATH and AJFM use the Ritt-Wu characteristic sets method, im- 
plemented by Wang (2001ab). The CharSets package, available in Maple (Wang, 
2002), is more versatile and powerful than our algorithm in Section IB~2l 

Lastly, Abbott et al. (2002) designed a Mathematica notebook with key func- 
tions for the computation of polynomial solutions in sn and cn . 
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Table 3: Comparison between codes PDESpecialSolutions.m and RATH. Test runs performed 
on a Dell Dimension 8200 PC, with 2.40 GHz Pentium 4 processor, 512 MB of RAM, with 
Mathcmatica v. 4.1 and Maple v. 7.0. The first 8 equations appear in Li and Liu (2002); the 
last 10 equations are listed in this paper. 

There are several symbolic tools for reducing and solving parameterized non- 
linear algebraic systems. Some are part of codes to simplify overdetermined ODE 
and PDE systems. For example, the Maple package Rif by Wittkopf and Reid 
(2003) allows for the computation of solution branches of nonlinear algebraic 
systems. The most powerful algebraic solvers use some flavor of the Grobner ba- 
sis algorithm. For up-to-date information on developments in this area we refer 
to Grabmeier et al. (2003). 

9. Discussion and Conclusions 

We presented several straightforward algorithms to compute special solutions of 
nonlinear PDEs, without using explicit integration. We designed the symbolic 
package PDESpecialSolutions.m to find solitary wave solutions of nonlinear 
PDEs involving tanh, sech, cn and sn functions. 

While the software reproduces the known (and also a few presumably new) so- 
lutions for many equations, there is no guarantee that the code will compute the 
complete solution set of all polynomial solutions involving the tanh and/or sech 
functions, especially when the PDEs have parameters. This is due to restrictions 
on the form of the solutions and the limitations of the algebraic solver. 

There is so much freedom in mixed tanh-sech solutions that the current code 
is limited to quadratic solutions. 

Furthermore, the nonlinear constraints which arise in solving the nonlinear 
algebraic system may be quintic or of higher degree, and therefore unsolvable in 
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analytic form. Also, since our software package is fully automated, it may not 
return the solutions in the simplest form. 

The example in Section 17131 illustrates this situation. By not solving quadratic 
or cubic equations explicitly the solutions (computed interactively with Mathe- 
matical can be presented in a more compact and readable form. 

In an attempt to avoid the explicit use of Mathematical Solve and Reduce 
functions, we considered various alternatives. For example, we used (i) variants of 
Grobner bases on the complete system, and (ii) combinatorics on the coefficients 
in the polynomial solutions (setting = or 7^ 0, for the admissible i and 
j). None of these alternatives paid off for systems with parameters. 

Often, the nonlinear solver returns constraints on the wave parameters Cj 
and the external parameters. In principle, one should verify whether or not such 
constraints affect the results of the previous steps in the algorithm. In particular, 
one should verify the consistency with the results from step 2 of the algorithms. 
We have not yet implemented this type of sophistication. 
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Appendix A: Using the software package 

We illustrate the use of the package PDESpecialSolutions .m on a PC. Users 
should have access to Mathematica v. 3.0 or higher. 

Put the package in a directory, say myDirectory, on drive C. Start a Mathe- 
matica notebook session and execute the commands: 

In[l] = SetDirectory ["c : WmyDirectory"] ; (* specify directory *) 

In[2] = <<PDESpecialSolutions .m (* read in package *) 

In [3] = PDESpecialSolutions [ 

{D[u[x,t] ,t] -alpha* (6*u[x, t] *D [u [x,t] ,x]+D[u[x,t] ,{x,3}])+ 
2*beta*v [x , t] *D [v [x , t] , x] == , 

D[v[x,t] ,t]+3*u[x,t] *D[v[x,t] ,x]+D[v[x,t] ,{x,3}] == 0}, 
{u[x,t] ,v[x,t] }, {x,t}, {alpha, beta}-, Form -> Sech, 
Verbose -> True, InputForm -> False, NumericTest -> True, 
SymbolicTest -> True, SolveAlgebraicSystem -> True 
(*, DegreeOfThePolynomial -> {m[l] -> 2, m[2] -> 1} *)] ; 

The package will compute the sech solutions (J3*7|) and (|3*%|) of the coupled KdV 
equation (jZ2j) - 

If the DegreeOfThePolynomial — > {m[l] — > 2, m[2] — > 1} were specified, the 
code would continue with this case only and compute (|3*7jl . 

If SolveAlgebraicSystem — > False, the algebraic system will be generated 
but not automatically solved. 

The format of PDESpecialSolutions is similar to the Mathematica function 
D Solve. The output is a list of sub-lists with solutions and constraints. The 
Backus-Naur Form of the function is: 

(Main Function) — > ~P'DESpecia.lSolutions[(E quations), (Functions), 

(Variables) , (Parameters), (Options)} 
(Options) — > Form — > (Form) | Verbose — > (Bool) \ 
InputForm — > (Bool) \ 

DegreeOfThePolynomial — > (List of Rules) \ 
SolveAlgebraicSystem — > (Bool) \ 
NumericTest — > (Bool) | SymbolicTest — > (Bool) 
(Form) — > Tanh | Sech | SechTanh | JacobiCN | JacobiSN 
(Bool) — > True | False 
(List of Rules) — > {m[l] — > Integer, m[2] — >• Integer, ...} 

The default value of Form is Tanh. The package PDESpecialSolutions .m has 

been tested on both UNIX work stations and PCs with Mathematica versions 
3.0, 4.0 and 4.1. A test set of over 50 PDEs and half a dozen ODEs was used. 



